###############################################################
## Liao and McDowell ISQ Analysis Replication R Code
## 9/25/2013
## R version 3.0.1 (2013-05-16) -- "Good Sport"
## NOTE: The following code ony works on Zelig version 3.5.5
##       as Zelig versions 4 and above have yet to implement 
##       rare events logit models combined with Amelia's 
##       multiple imputation.
###############################################################
# clear all objects in workspace
rm(list = ls())

# load libraries
library(foreign)
library(lubridate)
library(ggplot2)
library(RColorBrewer)
library(gridExtra)
library(zoo)
library(reshape2)
library(countrycode)
library(rworldmap)
library(rworldxtra)
library(plyr)
library(xtable)
library(reporttools)
library(Zelig)
library(MASS)
library(plm)
library(Amelia)
library(parallel)
library(texreg)

source("MI.table.r")
setwd("/Users/stevenliao/Dropbox/Research/replication/isq_rmb")
set.seed(1234)

# load data
data <- read.dta("yuan.dta")
data.settlement <- read.csv("yuan.settlement.csv", as.is=TRUE)
# 10 multiply imputed datasets in a list as a RData file
load("mi.rr.pt.RData")
# 10 multiply imputed datasets as separate Stata files
mi1 <- read.dta("mi.rr.pt1.dta")
mi2 <- read.dta("mi.rr.pt2.dta")
mi3 <- read.dta("mi.rr.pt3.dta")
mi4 <- read.dta("mi.rr.pt4.dta")
mi5 <- read.dta("mi.rr.pt5.dta")
mi6 <- read.dta("mi.rr.pt6.dta")
mi7 <- read.dta("mi.rr.pt7.dta")
mi8 <- read.dta("mi.rr.pt8.dta")
mi9 <- read.dta("mi.rr.pt9.dta")
mi10 <- read.dta("mi.rr.pt10.dta")

################################################
## Fig 1. RMB Settlement of Cross-Border Trade
################################################
data.settlement$trade_set_per1 <- data.settlement$trade_set_per*100
data.settlement.melt <- melt(data.settlement, measure.vars=c("trade_set_tot", "trade_set_per1"))
data.settlement.melt$variable <- factor(data.settlement.melt$variable, levels = c("trade_set_tot","trade_set_per1"), labels = c("RMB-based Trade (USD billions)", "% of China's Total Trade"))
p.ts3 <- ggplot(data.settlement.melt, aes(x = time, y=value)) + geom_bar(stat = "identity", width=.5) + facet_grid(variable~., scales="free_y", labeller = "label_value") + xlab("") + theme_bw()
# print graph
pdf("trade_settle.pdf", height = 6, width = 10)
p.ts3
dev.off()


#################################################################
## Fig 2. RMB Settlement of Outward and Inward Direct Investment
#################################################################
data.settlement.sub <- data.settlement[7:13,]
fdi.settle.long <- melt(data.settlement.sub, id = "time", measure = c("odi_set_tot", "fdi_set_tot"))
fdi.settle.long$time1 <- as.yearqtr(fdi.settle.long$time,format = "%YQ%q")
fdi.settle.long$time1 <- as.POSIXct(fdi.settle.long$time1)
# print graph
pdf("fdi_settle_bw.pdf", height = 4.5, width = 10)
ggplot(fdi.settle.long, aes(time1, value, colour = variable)) + geom_line() + xlab("") + ylab("RMB-based Direct Investment (USD billions)") + scale_colour_grey(name = "FDI Type", labels=c("ODI", "FDI")) + theme_bw()
dev.off()


#######################################################################################
## Fig 4. RMB Bilateral Swap Arrangement Countries by Swap Size (billions), 2009-2012
#######################################################################################
# subset by year
data.08 <- subset(data, year == 2008)
data.09 <- subset(data, year == 2009)
data.10 <- subset(data, year == 2010)
data.11 <- subset(data, year == 2011)

# join country data to map with rworldmap package
data.map.08 <- joinCountryData2Map(data.08, joinCode = "ISO3", nameJoinColumn = "iso3", verbose=TRUE)
data.map.09 <- joinCountryData2Map(data.09, joinCode = "ISO3", nameJoinColumn = "iso3", verbose=TRUE)
data.map.10 <- joinCountryData2Map(data.10, joinCode = "ISO3", nameJoinColumn = "iso3", verbose=TRUE)
data.map.11 <- joinCountryData2Map(data.11, joinCode = "ISO3", nameJoinColumn = "iso3", verbose=TRUE)

# define swap size breaks
swap.size.break <- c(0.1, 1, 50, 100, 200, 300, 400)

# generate black and white version of map
pdf("choropleth_swap_size_bw.pdf", height = 11, width = 8.5)#open pdf device
par(mfrow = c(4,1), mai=c(0,0,0.2,0.2))
# 2009
swap.map.09 <- mapCountryData(data.map.08, 
                              nameColumnToPlot = "size_lead", 
                              mapRegion = "world",
                              numCats = 6,
                              catMethod = swap.size.break,
                              missingCountryCol = "white",
                              borderCol = "lightgrey",
                              lwd = 0.2,
                              mapTitle = "2009",
                              addLegend = FALSE, 
                              colourPalette=brewer.pal(6, "Greys"))
do.call(addMapLegend, c(swap.map.09,
                        legendLabels = "all",
                        labelFontSize = 0.8,
                        legendMar = 4,
                        legendWidth = 0.5,
                        legendShrink = 0.9,
                        horizontal = FALSE))
# 2010
swap.map.10 <- mapCountryData(data.map.09, 
                              nameColumnToPlot = "size_lead", 
                              mapRegion = "world",
                              numCats = 6,
                              catMethod = swap.size.break,
                              missingCountryCol = "white",
                              borderCol = "lightgrey",
                              lwd = 0.2,
                              mapTitle = "2010",
                              addLegend = FALSE, 
                              colourPalette=brewer.pal(6, "Greys"))
do.call(addMapLegend, c(swap.map.10,
                        legendLabels = "all",
                        labelFontSize = 0.8,
                        legendMar = 4,
                        legendWidth = 0.5,
                        legendShrink = 0.9,
                        horizontal = FALSE))
# 2011
swap.map.11 <- mapCountryData(data.map.10, 
                              nameColumnToPlot = "size_lead", 
                              mapRegion = "world",
                              numCats = 6,
                              catMethod = swap.size.break,
                              missingCountryCol = "white",
                              borderCol = "lightgrey",
                              lwd = 0.2,
                              mapTitle = "2011",
                              addLegend = FALSE, 
                              colourPalette=brewer.pal(6, "Greys"))
do.call(addMapLegend, c(swap.map.11,
                        legendLabels = "all",
                        labelFontSize = 0.8,
                        legendMar = 4,
                        legendWidth = 0.5,
                        legendShrink = 0.9,
                        horizontal = FALSE))
# 2012
swap.map.12 <- mapCountryData(data.map.11, 
                              nameColumnToPlot = "size_lead", 
                              mapRegion = "world",
                              numCats = 6,
                              catMethod = swap.size.break,
                              missingCountryCol = "white",
                              borderCol = "lightgrey",
                              lwd = 0.2,
                              mapTitle = "2012",
                              addLegend = FALSE, 
                              colourPalette=brewer.pal(6, "Greys"))
do.call(addMapLegend, c(swap.map.12,
                        legendLabels = "all",
                        labelFontSize = 0.8,
                        legendMar = 4,
                        legendWidth = 0.5,
                        legendShrink = 0.9,
                        horizontal = FALSE))
dev.off()

#####################################
## Table 5. Descriptive Statistics
#####################################
# convert units
oil_prod_eia_mil <- data$oil_prod_eia/1000
coal_prod_eia_bil <- data$coal_prod_eia/1000000
oilexp_bil <- data$oilexp/1000000000
coalexpwt_bil <- data$coalexpwt/1000000000
copexpwt_mil <- data$copexpwt/1000000
irnexpwt_bil <- data$irnexpwt/1000000000
potexpwt_mil <- data$potexpwt/1000000
gdp2000_bil <- data$gdp2000/1000000000

# prepare for reporttools package
des.vars <- list(data$swap_lead,
                 data$ptradedep, data$chtradedep, 
                 data$pexpdep, data$pimpdep, data$chexpdep, data$chimpdep, 
                 data$pfdidep_book, data$chfdidep_book, 
                 data$tradeagree, data$bits_entry, 
                 data$scofull, 
                 oil_prod_eia_mil, coal_prod_eia_bil, 
                 oilexp_bil, coalexpwt_bil,
                 copexpwt_mil, irnexpwt_bil, potexpwt_mil,
                 gdp2000_bil, data$rgdpch, data$gdpgrow, data$dist_k)
des.nam <- c("Bilateral Swap Agreements (by end of 2012)",
             "Partner's Trade Dependence", "China's Trade Dependence",
             "Partner's Export Dependence", "Partner's Import Dependence", "China's Export Dependence", "China's Import Dependence", 
             "Partner's FDI Dependence", "China's FDI Dependence",
             "PTAs", "BITs",
             "SCO Member or Observer",
             "Oil Production (million barrels per day)","Coal Production (billion short tons)",
             "Oil Exports (billion USD)", "Coal Exports (billion kg)",
             "Copper Exports (million kg)", "Iron Exports (billion kg)", "Potash Exports (million kg)",
             "GDP (billion USD)", "GDP per capita (2005 international dollar)", "GDP Growth Rate", "Bilateral Distance (thousand km)")
des.tab.cont <-  tableContinuous(vars = des.vars,
                                 nams = des.nam,
                                 stats = c("mean", "min", "max", "n", "na"), 
                                 prec = 2, 
                                 col.tit.font = c("bf", "", "sf", "it", "rm"), 
                                 cap = "Descriptive Statistics", lab = "tab:des.tab",
                                 font.size = "normalsize", longtable = FALSE)


######################
## All Fitted models
######################
## Main models
# calculate tau for rare events logit
tau <- (sum(data$swap_lead))/(length(data$swap_lead)-sum(data$swap_lead))
# just key de facto variables with controls
relogit.robust.1 <- zelig(swap_lead ~ log(ptradedep)*log(chtradedep)
                          + log(pfdidep_book)*log(chfdidep_book)
                          + scofull 
                          + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                          + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                          model = "relogit", data = mi.more$imputations, 
                          tau = tau, case.control = "weighting", control = list(maxit = 500),
                          robust = TRUE)
relogit.robust.sum.1 <- summary(relogit.robust.1)

# just key de jure variables with controls
relogit.robust.2 <- zelig(swap_lead ~ tradeagree + bits_entry
                          + scofull 
                          + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                          + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                          model = "relogit", data = mi.more$imputations, 
                          tau = tau, case.control = "weighting", control = list(maxit = 500),
                          robust = TRUE)
relogit.robust.sum.2 <- summary(relogit.robust.2)

# full model (de facto + de jure)
relogit.robust.7 <- zelig(swap_lead ~ log(ptradedep)*log(chtradedep)
                          + log(pfdidep_book)*log(chfdidep_book)
                          + tradeagree + bits_entry
                          + scofull
                          + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                          + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                          model = "relogit", data = mi.more$imputations, 
                          tau = tau, case.control = "weighting", control = list(maxit = 500),
                          robust = TRUE)
relogit.robust.sum.7 <- summary(relogit.robust.7)

## Fitted robustness check models
# oil export instead of production
full.out.oilexp <- zelig(swap_lead ~ log(ptradedep)*log(chtradedep)
                         + log(pfdidep_book)*log(chfdidep_book)
                         + tradeagree + bits_entry
                         + scofull
                         + log(oilexp)
                         + log(coal_prod_eia+1)
                         + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                         model = "relogit", data = mi.more$imputations, 
                         tau = tau, case.control = "weighting", control = list(maxit = 500),
                         robust = TRUE)
full.sum.oilexp <- summary(full.out.oilexp)

# coal export instead of production
full.out.coalexp <- zelig(swap_lead ~ log(ptradedep)*log(chtradedep)
                          + log(pfdidep_book)*log(chfdidep_book)
                          + tradeagree + bits_entry
                          + scofull
                          + log(oil_prod_eia+1)
                          + log(coalexpwt)
                          + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                          model = "relogit", data = mi.more$imputations, 
                          tau = tau, case.control = "weighting", control = list(maxit = 500),
                          robust = TRUE)
full.sum.coalexp <- summary(full.out.coalexp)

# add copper
full.out.copper <- zelig(swap_lead ~ log(ptradedep)*log(chtradedep)
                         + log(pfdidep_book)*log(chfdidep_book)
                         + tradeagree + bits_entry
                         + scofull
                         + log(oil_prod_eia+1) + log(coal_prod_eia+1) + log(copexpwt)
                         + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                         model = "relogit", data = mi.more$imputations, 
                         tau = tau, case.control = "weighting", control = list(maxit = 500),
                         robust = TRUE)
full.sum.copper <- summary(full.out.copper)

# add iron
full.out.iron <- zelig(swap_lead ~ log(ptradedep)*log(chtradedep)
                       + log(pfdidep_book)*log(chfdidep_book)
                       + tradeagree + bits_entry
                       + scofull 
                       + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                       + log(irnexpwt)
                       + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                       model = "relogit", data = mi.more$imputations, 
                       tau = tau, case.control = "weighting", control = list(maxit = 500),
                       robust = TRUE)
full.sum.iron <- summary(full.out.iron)

# add potash
full.out.potash <- zelig(swap_lead ~ log(ptradedep)*log(chtradedep)
                         + log(pfdidep_book)*log(chfdidep_book)
                         + tradeagree + bits_entry
                         + scofull
                         + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                         + log(potexpwt)
                         + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                         model = "relogit", data = mi.more$imputations, 
                         tau = tau, case.control = "weighting", control = list(maxit = 500),
                         robust = TRUE)
full.sum.potash <- summary(full.out.potash) 

# Excluding Hong Kong
mi1.noHKG <- mi1[which(mi1$iso3!="HKG"), ]
mi2.noHKG <- mi2[which(mi2$iso3!="HKG"), ]
mi3.noHKG <- mi3[which(mi3$iso3!="HKG"), ]
mi4.noHKG <- mi4[which(mi4$iso3!="HKG"), ]
mi5.noHKG <- mi5[which(mi5$iso3!="HKG"), ]
mi6.noHKG <- mi6[which(mi6$iso3!="HKG"), ]
mi7.noHKG <- mi7[which(mi7$iso3!="HKG"), ]
mi8.noHKG <- mi8[which(mi8$iso3!="HKG"), ]
mi9.noHKG <- mi9[which(mi9$iso3!="HKG"), ]
mi10.noHKG <- mi10[which(mi10$iso3!="HKG"), ]
# full model without HKG
full.out.noHKG <- zelig(swap_lead ~ log(ptradedep)*log(chtradedep)
                        + log(pfdidep_book)*log(chfdidep_book)
                        + tradeagree + bits_entry
                        + scofull
                        + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                        + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                        model = "relogit", data = mi(mi1.noHKG, mi2.noHKG, mi3.noHKG, mi4.noHKG, mi5.noHKG, mi6.noHKG, mi7.noHKG, mi8.noHKG, mi9.noHKG, mi10.noHKG), 
                        tau = tau, case.control = "weighting", control = list(maxit = 500),
                        robust = TRUE)
full.sum.noHKG <- summary(full.out.noHKG)

# Excluding Macao
mi1.noMAC <- mi1[which(mi1$iso3!="MAC"), ]
mi2.noMAC <- mi2[which(mi2$iso3!="MAC"), ]
mi3.noMAC <- mi3[which(mi3$iso3!="MAC"), ]
mi4.noMAC <- mi4[which(mi4$iso3!="MAC"), ]
mi5.noMAC <- mi5[which(mi5$iso3!="MAC"), ]
mi6.noMAC <- mi6[which(mi6$iso3!="MAC"), ]
mi7.noMAC <- mi7[which(mi7$iso3!="MAC"), ]
mi8.noMAC <- mi8[which(mi8$iso3!="MAC"), ]
mi9.noMAC <- mi9[which(mi9$iso3!="MAC"), ]
mi10.noMAC <- mi10[which(mi10$iso3!="MAC"), ]
# full model without MAC
full.out.noMAC <- zelig(swap_lead ~ log(ptradedep)*log(chtradedep)
                        + log(pfdidep_book)*log(chfdidep_book)
                        + tradeagree + bits_entry
                        + scofull
                        + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                        + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                        model = "relogit", data = mi(mi1.noMAC, mi2.noMAC, mi3.noMAC, mi5.noMAC, mi5.noMAC, mi6.noMAC, mi7.noMAC, mi8.noMAC, mi9.noMAC, mi10.noMAC), 
                        tau = tau, case.control = "weighting", control = list(maxit = 500),
                        robust = TRUE)
full.sum.noMAC <- summary(full.out.noMAC)

# export interdependence
full.out.exdep <- zelig(swap_lead ~ log(pexpdep)*log(chexpdep)
                        + log(pimpdep) + log(chimpdep)
                        + log(pfdidep_book)*log(chfdidep_book)
                        + tradeagree + bits_entry
                        + scofull
                        + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                        + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                        model = "relogit", data = mi.more$imputations, 
                        tau = tau, case.control = "weighting", control = list(maxit = 500),
                        robust = TRUE)
full.sum.exdep <- summary(full.out.exdep)

# import interdependence
full.out.imdep <- zelig(swap_lead ~ log(pexpdep) + log(chexpdep)
                        + log(pimpdep)*log(chimpdep)
                        + log(pfdidep_book)*log(chfdidep_book)
                        + tradeagree + bits_entry
                        + scofull
                        + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                        + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                        model = "relogit", data = mi.more$imputations, 
                        tau = tau, case.control = "weighting", control = list(maxit = 500),
                        robust = TRUE)
full.sum.imdep <- summary(full.out.imdep)

# pimpdep*chexpdep
full.out.pimp.chexp <- zelig(swap_lead ~ log(pexpdep) + log(chimpdep)
                             + log(pimpdep)*log(chexpdep)
                             + log(pfdidep_book)*log(chfdidep_book)
                             + tradeagree + bits_entry
                             + scofull
                             + log(oil_prod_eia+1) + log(coal_prod_eia+1)
                             + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                             model = "relogit", data = mi.more$imputations, 
                             tau = tau, case.control = "weighting", control = list(maxit = 500),
                             robust = TRUE)
full.sum.pimp.chexp <- summary(full.out.pimp.chexp)

# pexpdep*chimpdep
full.out.pexp.chimp <- zelig(swap_lead ~ log(pimpdep) + log(chexpdep)
                             + log(pexpdep)*log(chimpdep)
                             + log(pfdidep_book)*log(chfdidep_book)
                             + tradeagree + bits_entry
                             + scofull
                             + log(oil_prod_eia+1) + log(coal_prod_eia+1) 
                             + log(gdp2000) + log(rgdpch) + gdpgrow  + dist_k,
                             model = "relogit", data = mi.more$imputations, 
                             tau = tau, case.control = "weighting", control = list(maxit = 500),
                             robust = TRUE)
full.sum.pexp.chimp <- summary(full.out.pexp.chimp)

# dynamic probit
d.prob.7 <- zelig(swap_lead ~ I(swap)
                  + log(ptradedep)*log(chtradedep) + I(swap*log(ptradedep))*I(swap*log(chtradedep))
                  + log(pfdidep_book)*log(chfdidep_book) + I(swap*log(pfdidep_book))*I(swap*log(chfdidep_book))
                  + tradeagree + I(swap*tradeagree) 
                  + bits_entry + I(swap*bits_entry)
                  + scofull + I(swap*scofull)
                  + log(oil_prod_eia+1) + I(swap*log(oil_prod_eia+1))
                  + log(coal_prod_eia+1) + I(swap*log(coal_prod_eia+1))
                  + log(gdp2000) + I(swap*log(gdp2000))
                  + log(rgdpch) + I(swap*log(rgdpch))
                  + gdpgrow + I(swap*gdpgrow)
                  + dist_k + I(swap*dist_k),
                  model = "probit", data = mi(mi1,mi2,mi3,mi4,mi5,mi6,mi7,mi8,mi9,mi10), control = list(maxit = 500),
                  robust = TRUE)
d.prob.sum.7 <- summary(d.prob.7)


##############################
## Fig 6. Coefficients Plot
##############################
# de facto model
coef.names.1 <- rownames(coef(relogit.robust.sum.1))
coef.values.1 <- coef(relogit.robust.sum.1)[,1]
se.values.1 <- coef(relogit.robust.sum.1)[,2]
# de jure model
coef.names.2 <- rownames(coef(relogit.robust.sum.2))
coef.values.2 <- coef(relogit.robust.sum.2)[,1]
se.values.2 <- coef(relogit.robust.sum.2)[,2]
# full model
coef.names.7 <- rownames(coef(relogit.robust.sum.7))
coef.values.7 <- coef(relogit.robust.sum.7)[,1]
se.values.7 <- coef(relogit.robust.sum.7)[,2]

# append NAs to have same vector length as full model
coef.append.1 <- append(coef.values.1, c(NA,NA), after=5)
se.append.1 <- append(se.values.1, c(NA,NA), after=5)
coef.append.2 <- append(coef.values.2, c(NA,NA,NA,NA), after=1)
se.append.2 <- append(se.values.2, c(NA,NA,NA,NA), after=1)
coef.append.2 <- append(coef.append.2, c(NA,NA), after=14)
se.append.2 <- append(se.append.2, c(NA,NA), after=14)

# remove intercepts
coef.vec.1 <- coef.append.1[-1]
se.vec.1 <- se.append.1[-1]
coef.vec.2 <- coef.append.2[-1]
se.vec.2 <- se.append.2[-1]
coef.vec.7 <- coef.values.7[-1]
se.vec.7 <- se.values.7[-1]

# list appropriate variable names
var.names <- c("Log(Partner's Trade Dependence)", "Log(China's Trade Dependence)",
               "Log(Partner's FDI Dependence)","Log(China's FDI Dependence)",  
               "PTAs","BITs", 
               "SCO", 
               "Log(Oil Production)", "Log(Coal Production)",
               "Log(GDP)", "Log(GDP per capita)", "GDP Growth", "Distance (Thousand Km)",
               "Trade Interdependence", "FDI Interdependence")

# prepare parameters for plot
y.axis <- c(length(coef.vec.7):1)#create indicator for y.axis, descending so that R orders vars from top to bottom on y-axis
adjust <- .15#create object that we will use to adjust points and lines up and down to distinguish between models
pdf("coefplot.pdf", height = 7, width = 7)#open pdf device
par(mfrow = c(1,1))#reset graphical window
par(mar = c(2, 12, .5, 0))#set margins for plot, leaving lots of room on left-margin (2nd number in margin command) for variable names
plot(coef.vec.1, y.axis + adjust, type = "p", axes = F, xlab = "", ylab = "", pch = 17, cex = .6, col = "black",#plot coefficients as points, turning off axes and labels. 
     xlim = c(-1,5), ylim = c(min(y.axis), max(y.axis)), xaxs = "r", main = "") #set limits of x-axis so that they include mins and maxs of 
#coefficients + 95% confidence intervals and plot is symmetric; use "internal axes", and leave plot title empty
segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis+adjust, coef.vec.1+qnorm(.975)*se.vec.1, y.axis+adjust, lwd =  1.5, col="black")#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
axis(1, at = seq(-1,5,by=1), labels =  seq(-1,5,by=1), tick = T,#draw x-axis and labels with tick marks
     cex.axis = .8, mgp = c(2,.5,0))#reduce label size, moves labels closer to tick marks
axis(2, at = y.axis, label = var.names, las = 1, tick = T, 
     cex.axis = .8) #draw y-axis with tick marks, make labels perpendicular to axis and closer to axis
segments(0,0,0,17,lty=2) # draw dotted line through 0

# add de jure model
#because we are using white points and do want the lines to go "through" points rather than over them 
#we draw lines first and the overlay points
segments(coef.vec.2-qnorm(.975)*se.vec.2, y.axis, coef.vec.2+qnorm(.975)*se.vec.2, y.axis, lwd =  1.3, col = "black")#draw lines connecting 95% confidence intervals
points(coef.vec.2, y.axis, pch = 24, cex = 0.6, col = "black", bg = "white") #add point estimates for 2nd model; pch = 24 uses for overlay points, and "white" for white color

# add full model
segments(coef.vec.7-qnorm(.975)*se.vec.7, y.axis-adjust, coef.vec.7+qnorm(.975)*se.vec.7, y.axis-adjust, lwd =  1.3, col = "black")#draw lines connecting 95% confidence intervals
points(coef.vec.7, y.axis-adjust, pch = 19, cex = 0.6, col = "black" ) #add point estimates for 3nd model; pch = 19 uses for overlay points, and "white" for white color

# add legend (manually) to identify which dots denote model 1, 2, and 3
points(2.6, 3, pch = 17, cex  = 0.6, col = "black")
text(2.8, 3, "De facto Model", pos = 4, cex  = 0.8)
points(2.6, 2.5, pch = 24, cex  = 0.6, col = "black")
text(2.8, 2.5, "De jure Model", pos = 4, cex  = 0.8)
points(2.6, 2, pch = 19, cex  = 0.6, col = "black")
text(2.8, 2, "Full Model", pos = 4, cex  = 0.8)
segments(2.6, 1.5, 2.8, 1.5, lwd = 1, col = "black")
text(2.8, 1.5, "95% C.I.", pos = 4, cex  = 0.8)
dev.off() #turn off pdf device; graph is created in working directory


#############################################################
## Table 6. Parameter Estimates and Robust Standard Errors
#############################################################
## prepare vectors for texreg package 
# de facto model
coef.names.1 <- rownames(coef(relogit.robust.sum.1))
coef.values.1 <- coef(relogit.robust.sum.1)[,1]
se.values.1 <- coef(relogit.robust.sum.1)[,2]
p.val.1 <- coef(relogit.robust.sum.1)[,4]
aic.1 <- (relogit.robust.1[[1]][[11]]+relogit.robust.1[[2]][[11]]+relogit.robust.1[[3]][[11]]+relogit.robust.1[[4]][[11]]+relogit.robust.1[[5]][[11]]+relogit.robust.1[[6]][[11]]+relogit.robust.1[[7]][[11]]+relogit.robust.1[[8]][[11]]+relogit.robust.1[[9]][[11]]+relogit.robust.1[[10]][[11]])/10

# de jure model
coef.names.2 <- rownames(coef(relogit.robust.sum.2))
coef.values.2 <- coef(relogit.robust.sum.2)[,1]
se.values.2 <- coef(relogit.robust.sum.2)[,2]
p.val.2 <- coef(relogit.robust.sum.2)[,4]
aic.2 <- (relogit.robust.2[[1]][[11]]+relogit.robust.2[[2]][[11]]+relogit.robust.2[[3]][[11]]+relogit.robust.2[[4]][[11]]+relogit.robust.2[[5]][[11]]+relogit.robust.2[[6]][[11]]+relogit.robust.2[[7]][[11]]+relogit.robust.2[[8]][[11]]+relogit.robust.2[[9]][[11]]+relogit.robust.2[[10]][[11]])/10

# full model
coef.names.7 <- rownames(coef(relogit.robust.sum.7))
coef.values.7 <- coef(relogit.robust.sum.7)[,1]
se.values.7 <- coef(relogit.robust.sum.7)[,2]
p.val.7 <- coef(relogit.robust.sum.7)[,4]
aic.7 <- (relogit.robust.7[[1]][[11]]+relogit.robust.7[[2]][[11]]+relogit.robust.7[[3]][[11]]+relogit.robust.7[[4]][[11]]+relogit.robust.7[[5]][[11]]+relogit.robust.7[[6]][[11]]+relogit.robust.7[[7]][[11]]+relogit.robust.7[[8]][[11]]+relogit.robust.7[[9]][[11]]+relogit.robust.7[[10]][[11]])/10

n <- length(mi1$swap)

defacto <- MI(names=coef.names.1, 
              coef=coef.values.1, 
              se=se.values.1, 
              pval=p.val.1, 
              n=n,
              aic=aic.1)
dejure <- MI(names=coef.names.2, 
             coef=coef.values.2, 
             se=se.values.2, 
             pval=p.val.2, 
             n=n,
             aic=aic.2)
full <- MI(names=coef.names.7, 
           coef=coef.values.7, 
           se=se.values.7, 
           pval=p.val.7, 
           n=n,
           aic=aic.7)

## make table with texreg
texreg(list(defacto,dejure,full), custom.model.names=c("De facto", "De jure", "Full Relogit"), 
       stars = c(0.001, 0.01, 0.05, 0.1), symbol="\\wedge", 
       custom.coef.names = c("Intercept", "log(P. Trade Dependence)", "log(CN Trade Dependence)",
                             "log(P. FDI Dep.)", "log(CN FDI Dep.)", "SCO","log(Oil Production)",
                             "log(Coal Production)", "log(GDP)", "log(GDP per capita)", "GDP Growth Rate",
                             "Bilateral Distance", "log(P. Trade Dep.):log(CN Trade Dep.)", 
                             "log(P. FDI Dep.):log(CN FDI Dep.)","PTAs","BITs"), 
       reorder.coef = c(1, 2, 3, 13, 4, 5, 14, 15, 16, 6, 7, 8, 9, 10, 11, 12), scriptsize = TRUE, single.row = TRUE, 
       float.pos = "!ht", caption = "Parameter Estimates and Robust Standard Errors", label = "Estimates", 
       use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE)

#######################################
## Simulated Predicted Probabilities
#######################################
# set trade ranges and values for mean -/+ 1 s.d.
ptradedep.range <- with(data, seq(min(ptradedep, na.rm = TRUE), max(ptradedep, na.rm = TRUE), length.out = 1000))
chtradedep.msd <- ifelse(mean(data$chtradedep, na.rm = TRUE) - sd(data$chtradedep, na.rm = TRUE) < 0,
                         min(data$chtradedep, na.rm=TRUE),
                         mean(data$chtradedep, na.rm = TRUE) - sd(data$chtradedep, na.rm = TRUE) < 0)
chtradedep.psd <- mean(data$chtradedep, na.rm = TRUE) + sd(data$chtradedep, na.rm = TRUE)
chtradedep.p2sd <- mean(data$chtradedep, na.rm = TRUE) + sd(data$chtradedep, na.rm = TRUE)*2
chtradedep.range <- with(data, seq(min(chtradedep, na.rm = TRUE), max(chtradedep, na.rm = TRUE), length.out = 1000))
ptradedep.msd <- ifelse(mean(data$ptradedep, na.rm = TRUE) - sd(data$ptradedep, na.rm = TRUE) < 0,
                        min(data$ptradedep, na.rm = TRUE),
                        mean(data$ptradedep, na.rm = TRUE) - sd(data$ptradedep, na.rm = TRUE))
ptradedep.psd <- mean(data$ptradedep, na.rm = TRUE) + sd(data$ptradedep, na.rm = TRUE)
ptradedep.p2sd <- mean(data$ptradedep, na.rm = TRUE) + sd(data$ptradedep, na.rm = TRUE)*2

# simulate predicted probabilities with ptradedep as x axis and chtradedep as first difference
chtradedep.low <- setx(relogit.robust.7, chtradedep=chtradedep.msd, ptradedep = ptradedep.range, scofull=0, tradeagree=0, bits_entry=0)
chtradedep.high <- setx(relogit.robust.7, chtradedep=chtradedep.psd, ptradedep = ptradedep.range, scofull=0, tradeagree=0,bits_entry=0)
s.out <-sim(relogit.robust.7, x = chtradedep.low, x1 = chtradedep.high, num = 1000)

# simulate predicted probabilities with chtradedep as x axis and ptradedep as first difference
ptradedep.low <- setx(relogit.robust.7, chtradedep=chtradedep.range, ptradedep = ptradedep.msd, scofull=0, tradeagree=0, bits_entry=0)
ptradedep.high <- setx(relogit.robust.7, chtradedep=chtradedep.range, ptradedep = ptradedep.psd, scofull=0, tradeagree=0, bits_entry=0)
s.out.2 <-sim(relogit.robust.7, x = ptradedep.low, x1 = ptradedep.high, num = 1000)

# get fd for ptradedep when chtradedep is high
ptradedep.low.1 <- setx(relogit.robust.7, chtradedep=chtradedep.psd, ptradedep = ptradedep.msd, scofull=0, tradeagree=0, bits_entry=0)
ptradedep.high.1 <- setx(relogit.robust.7, chtradedep=chtradedep.psd, ptradedep = ptradedep.psd, scofull=0, tradeagree=0, bits_entry=0)
s.out.3 <-sim(relogit.robust.7, x = ptradedep.low.1, x1 = ptradedep.high.1, num = 1000)
s.out.ptradedep.fd.stats <- summary(s.out.3)$qi.stats$fd
ptradedep.fd.plot <- s.out.ptradedep.fd.stats[1]

# get fd for chtradedep when ptradedep is high
chtradedep.low.1 <- setx(relogit.robust.7, chtradedep=chtradedep.msd, ptradedep = ptradedep.psd, scofull=0, tradeagree=0, bits_entry=0)
chtradedep.high.1 <- setx(relogit.robust.7, chtradedep=chtradedep.psd, ptradedep = ptradedep.psd, scofull=0, tradeagree=0, bits_entry=0)
s.out.4 <-sim(relogit.robust.7, x = chtradedep.low.1, x1 = chtradedep.high.1, num = 1000)
s.out.chtradedep.fd.stats <- summary(s.out.4)$qi.stats$fd
chtradedep.fd.plot <- s.out.chtradedep.fd.stats[1]

# set ranges and values for mean -/+ 1 s.d.
pfdidep_book.range <- with(data, seq(min(pfdidep_book, na.rm = TRUE), max(pfdidep_book, na.rm = TRUE), length.out = 1000))
chfdidep_book.msd <- ifelse(mean(data$chfdidep_book, na.rm = TRUE) - sd(data$chfdidep_book, na.rm = TRUE) < 0,
                            min(data$chfdidep_book, na.rm=TRUE),
                            mean(data$chfdidep_book, na.rm = TRUE) - sd(data$chfdidep_book, na.rm = TRUE))
chfdidep_book.psd <- mean(data$chfdidep_book, na.rm = TRUE) + sd(data$chfdidep_book, na.rm = TRUE)
chfdidep_book.p2sd <- mean(data$chfdidep_book, na.rm = TRUE) + sd(data$chfdidep_book, na.rm = TRUE)*2
chfdidep_book.range <- with(data, seq(min(chfdidep_book, na.rm = TRUE), max(chfdidep_book, na.rm = TRUE), length.out = 1000))
pfdidep_book.msd <- ifelse(mean(data$pfdidep_book, na.rm = TRUE) - sd(data$pfdidep_book, na.rm = TRUE) < 0,
                           min(data$pfdidep_book, na.rm=TRUE), 
                           mean(data$pfdidep_book, na.rm = TRUE) - sd(data$pfdidep_book, na.rm = TRUE))
pfdidep_book.psd <- mean(data$pfdidep_book, na.rm = TRUE) + sd(data$pfdidep_book, na.rm = TRUE)
pfdidep_book.p2sd <- mean(data$pfdidep_book, na.rm = TRUE) + sd(data$pfdidep_book, na.rm = TRUE)*2

# simulate predicted probabilities with pfdidep_book as x axis and chfdidep_book as first difference
chfdidep_book.low <- setx(relogit.robust.7, chfdidep_book=chfdidep_book.msd, pfdidep_book = pfdidep_book.range, scofull=0, tradeagree=0, bits_entry=0)
chfdidep_book.high <- setx(relogit.robust.7, chfdidep_book=chfdidep_book.psd, pfdidep_book = pfdidep_book.range, scofull=0, tradeagree=0, bits_entry=0)
s.out.fdi.px <-sim(relogit.robust.7, x = chfdidep_book.low, x1 = chfdidep_book.high, num = 1000)

# simulate predicted probabilities with chfdidep_book as x axis and pfdidep_book as first difference
pfdidep_book.low <- setx(relogit.robust.7, pfdidep_book=pfdidep_book.msd, chfdidep_book = chfdidep_book.range, scofull=0, tradeagree=0, bits_entry=0)
pfdidep_book.high <- setx(relogit.robust.7, pfdidep_book=pfdidep_book.psd, chfdidep_book = chfdidep_book.range, scofull=0, tradeagree=0,bits_entry=0)
s.out.fdi.cx <-sim(relogit.robust.7, x = pfdidep_book.low, x1 = pfdidep_book.high, num = 1000)

# get fd for pfdidep_book when chfdidep_book is high
pfdidep_book.low.1 <- setx(relogit.robust.7, chfdidep_book = chfdidep_book.psd, pfdidep_book = pfdidep_book.msd, scofull=0, tradeagree=0, bits_entry=0)
pfdidep_book.high.1 <- setx(relogit.robust.7, chfdidep_book = chfdidep_book.psd, pfdidep_book = pfdidep_book.psd, scofull=0, tradeagree=0, bits_entry=0)
s.out.fdi.px.1 <-sim(relogit.robust.7, x = pfdidep_book.low.1, x1 = pfdidep_book.high.1, num = 1000)
s.out.pfdidep_book.fd.stats <- summary(s.out.fdi.px.1)$qi.stats$fd
pfdidep_book.fd.plot <- s.out.pfdidep_book.fd.stats[1]

# get fd for chfdidep_book when pfdidep_book is high
chfdidep_book.low.1 <- setx(relogit.robust.7, pfdidep_book = pfdidep_book.psd, chfdidep_book = chfdidep_book.msd, scofull=0, tradeagree=0, bits_entry=0)
chfdidep_book.high.1 <- setx(relogit.robust.7, pfdidep_book = pfdidep_book.psd, chfdidep_book = chfdidep_book.psd, scofull=0, tradeagree=0, bits_entry=0)
s.out.fdi.cx.1 <-sim(relogit.robust.7, x = chfdidep_book.low.1, x1 = chfdidep_book.high.1, num = 1000)
s.out.chfdidep_book.fd.stats <- summary(s.out.fdi.cx.1)$qi.stats$fd
chfdidep_book.fd.plot <- s.out.chfdidep_book.fd.stats[1]

## other simulated first differences
# simulate first difference for PTAs
pta.no <- setx(relogit.robust.7, tradeagree=0, bits_entry=0, scofull=0)
pta.yes <- setx(relogit.robust.7, tradeagree=1, bits_entry=0, scofull=0)
s.out.pta <- sim(relogit.robust.7, x = pta.no, x1 = pta.yes, num = 1000)
s.out.pta.fd.stats <- summary(s.out.pta)$qi.stats$fd
pta.fd.plot <- s.out.pta.fd.stats[1]
pta.fd <- pta.fd.plot*100

# simulate first difference for BITs
bits.no <- setx(relogit.robust.7, tradeagree=0, bits_entry=0, scofull=0)
bits.yes <- setx(relogit.robust.7, tradeagree=0, bits_entry=1, scofull=0)
s.out.bits <- sim(relogit.robust.7, x = bits.no, x1 = bits.yes, num = 1000)
s.out.bits.fd.stats <- summary(s.out.bits)$qi.stats$fd
bits.fd.plot <- s.out.bits.fd.stats[1]
bits.fd <- bits.fd.plot*100

# simulate first difference for SCO
sco.no <- setx(relogit.robust.7, tradeagree=0, bits_entry=0, scofull=0)
sco.yes <- setx(relogit.robust.7, tradeagree=0, bits_entry=0, scofull=1)
s.out.sco <- sim(relogit.robust.7, x = sco.no, x1 = sco.yes, num = 1000)
s.out.sco.fd.stats <- summary(s.out.sco)$qi.stats$fd
sco.fd.plot <- s.out.sco.fd.stats[1]
sco.fd <- sco.fd.plot*100

# simulate first difference for rgdpch
rgdpch.msd <- ifelse(mean(data$rgdpch, na.rm = TRUE) - sd(data$rgdpch, na.rm = TRUE) < 0,
                     min(data$rgdpch, na.rm = TRUE),
                     mean(data$rgdpch, na.rm = TRUE) - sd(data$rgdpch, na.rm = TRUE))
rgdpch.psd <- mean(data$rgdpch, na.rm = TRUE) + sd(data$rgdpch, na.rm = TRUE)
rgdpch.low <- setx(relogit.robust.7, rgdpch = rgdpch.msd, tradeagree=0, bits_entry=0, scofull=0)
rgdpch.high <- setx(relogit.robust.7, rgdpch = rgdpch.psd, tradeagree=0, bits_entry=0, scofull=0)
s.out.rgdpch <- sim(relogit.robust.7, x = rgdpch.low, x1 = rgdpch.high, num = 1000)
s.out.rgdpch.fd.stats <- summary(s.out.rgdpch)$qi.stats$fd
rgdpch.fd.plot <- s.out.rgdpch.fd.stats[1]
rgdpch.fd <- rgdpch.fd.plot*100

# simulate first difference for gdpgrow
gdpgrow.msd <- ifelse(mean(data$gdpgrow, na.rm = TRUE) - sd(data$gdpgrow, na.rm = TRUE) < 0,
                      min(data$gdpgrow, na.rm = TRUE),
                      mean(data$gdpgrow, na.rm = TRUE) - sd(data$gdpgrow, na.rm = TRUE))
gdpgrow.psd <- mean(data$gdpgrow, na.rm = TRUE) + sd(data$gdpgrow, na.rm = TRUE)
gdpgrow.low <- setx(relogit.robust.7, gdpgrow = gdpgrow.msd, tradeagree=0, bits_entry=0, scofull=0)
gdpgrow.high <- setx(relogit.robust.7, gdpgrow = gdpgrow.psd, tradeagree=0, bits_entry=0, scofull=0)
s.out.gdpgrow <- sim(relogit.robust.7, x = gdpgrow.low, x1 = gdpgrow.high, num = 1000)
s.out.gdpgrow.fd.stats <- summary(s.out.gdpgrow)$qi.stats$fd
gdpgrow.fd.plot <- s.out.gdpgrow.fd.stats[1]
gdpgrow.fd <- gdpgrow.fd.plot*100

# simulate first difference for oil_prod_eia
oil_prod_eia.msd <- ifelse(mean(data$oil_prod_eia, na.rm = TRUE) - sd(data$oil_prod_eia, na.rm = TRUE) < 0,
                           min(data$oil_prod_eia, na.rm = TRUE),
                           mean(data$oil_prod_eia, na.rm = TRUE) - sd(data$oil_prod_eia, na.rm = TRUE))
oil_prod_eia.psd <- mean(data$oil_prod_eia, na.rm = TRUE) + sd(data$oil_prod_eia, na.rm = TRUE)
oil_prod_eia.low <- setx(relogit.robust.7, oil_prod_eia = oil_prod_eia.msd, tradeagree=0, bits_entry=0, scofull=0)
oil_prod_eia.high <- setx(relogit.robust.7, oil_prod_eia = oil_prod_eia.psd, tradeagree=0, bits_entry=0, scofull=0)
s.out.oil_prod_eia <- sim(relogit.robust.7, x = oil_prod_eia.low, x1 = oil_prod_eia.high, num = 1000)
s.out.oil_prod_eia.fd.stats <- summary(s.out.oil_prod_eia)$qi.stats$fd
oil_prod_eia.fd.plot <- s.out.oil_prod_eia.fd.stats[1]
oil_prod_eia.fd <- oil_prod_eia.fd.plot*100


# simulate first difference for coal_prod_eia
coal_prod_eia.msd <- ifelse(mean(data$coal_prod_eia, na.rm = TRUE) - sd(data$coal_prod_eia, na.rm = TRUE) < 0,
                            min(data$coal_prod_eia, na.rm = TRUE),
                            mean(data$coal_prod_eia, na.rm = TRUE) - sd(data$coal_prod_eia, na.rm = TRUE) < 0)
coal_prod_eia.psd <- mean(data$coal_prod_eia, na.rm = TRUE) + sd(data$coal_prod_eia, na.rm = TRUE)
coal_prod_eia.low <- setx(relogit.robust.7, coal_prod_eia = coal_prod_eia.msd, tradeagree=0, bits_entry=0, scofull=0)
coal_prod_eia.high <- setx(relogit.robust.7, coal_prod_eia = coal_prod_eia.psd, tradeagree=0, bits_entry=0, scofull=0)
s.out.coal_prod_eia <- sim(relogit.robust.7, x = coal_prod_eia.low, x1 = coal_prod_eia.high, num = 1000)
s.out.coal_prod_eia.fd.stats <- summary(s.out.coal_prod_eia)$qi.stats$fd
coal_prod_eia.fd.plot <- s.out.coal_prod_eia.fd.stats[1]
coal_prod_eia.fd <- coal_prod_eia.fd.plot*100

# simulate first difference for gdp2000
gdp2000.msd <- ifelse(mean(data$gdp2000, na.rm = TRUE) - sd(data$gdp2000, na.rm = TRUE) < 0, 
                      min(data$gdp2000, na.rm = TRUE), 
                      mean(data$gdp2000, na.rm = TRUE) - sd(data$gdp2000, na.rm = TRUE))
gdp2000.psd <- mean(data$gdp2000, na.rm = TRUE) + sd(data$gdp2000, na.rm = TRUE)
gdp2000.low <- setx(relogit.robust.7, gdp2000 = gdp2000.msd, tradeagree=0, bits_entry=0, scofull=0)
gdp2000.high <- setx(relogit.robust.7, gdp2000 = gdp2000.psd, tradeagree=0, bits_entry=0, scofull=0)
s.out.gdp2000 <- sim(relogit.robust.7, x = gdp2000.low, x1 = gdp2000.high, num = 1000)
s.out.gdp2000.fd.stats <- summary(s.out.gdp2000)$qi.stats$fd
gdp2000.fd.plot <- s.out.gdp2000.fd.stats[1]
gdp2000.fd <- gdp2000.fd.plot*100

# simulate first difference for dist_k
dist_k.msd <- ifelse(mean(data$dist_k, na.rm = TRUE) - sd(data$dist_k, na.rm = TRUE) < 0, 
                     min(data$dist_k, na.rm = TRUE), 
                     mean(data$dist_k, na.rm = TRUE) - sd(data$dist_k, na.rm = TRUE))
dist_k.psd <- mean(data$dist_k, na.rm = TRUE) + sd(data$dist_k, na.rm = TRUE)
dist_k.low <- setx(relogit.robust.7, dist_k = dist_k.msd, tradeagree=0, bits_entry=0, scofull=0)
dist_k.high <- setx(relogit.robust.7, dist_k = dist_k.psd, tradeagree=0, bits_entry=0, scofull=0)
s.out.dist_k <- sim(relogit.robust.7, x = dist_k.low, x1 = dist_k.high, num = 1000)
s.out.dist_k.fd.stats <- summary(s.out.dist_k)$qi.stats$fd
dist_k.fd.plot <- s.out.dist_k.fd.stats[1]
dist_k.fd <- dist_k.fd.plot*100


########################################################################################################
## Fig 7. Partner’s Trade Dependence on China and Simulated Predicted Probabilities of BSA Occurrences
## ptradedep as x axis, chtradedep as first difference
########################################################################################################
# make a tall data.frame from the Zelig simulation object
toReshape_px <- data.frame(x = s.out[[2]][, "log(ptradedep)"], t(s.out[[6]]$ev))
longEV_px <- melt(toReshape_px, id.vars = "x")
toReshape2_px <- data.frame(x = s.out[[2]][, "log(ptradedep)"], t(s.out[[6]]$fd) + t(s.out[[6]]$ev))
longFD_px <- melt(toReshape2_px, id.vars = "x")
longEV_px$Setting <- "Low Chinese Trade Dependence on Partner"
longFD_px$Setting <- "High Chinese Trade Dependence on Partner"
longEV_px <- data.frame(rbind(longEV_px, longFD_px))
longEV_px$x_percent <- exp(longEV_px$x)*100

# plot
zp.px.facet <- ggplot(longEV_px, aes(x = x_percent, y = value, group = variable)) +
  geom_line(alpha = I(1/sqrt(nrow(s.out[[6]]$ev)))) +
  #geom_point(shape = 20, color = "gray30", alpha = I(0.05)) +
  #stat_smooth(method = "gam", se = FALSE, formula = y ~ poly(x, 3)) +
  facet_grid(Setting~., scales="free_y", labeller = "label_value") +
  scale_x_log10("Partner's Trade Dependence on China (%), scale = log_10",
                breaks=c(0.1, 1, 10, 15), 
                labels=c("0.1", "1", "10", "15"), 
                expand = c(0, 0)) +
  scale_y_continuous("Predicted Probability of BSAs", limits = c(0, 1), expand = c(0, 0.01)) +
  theme_bw()
pdf(file="predictedprob_facet_bw.pdf", width=7, height=7)
print(zp.px.facet)
dev.off()

########################################################################################################
## Fig 8. Partner’s FDI Dependence on China and Simulated Predicted Probabilities of BSA Occurrences
## pfdidep as x axis, chfdidep as first difference
########################################################################################################
# Make a tall data.frame from the Zelig simulation object
toReshape.fdi.px <- data.frame(x = s.out.fdi.px[[2]][, "log(pfdidep_book)"], t(s.out.fdi.px[[6]]$ev))
longEV.fdi.px <- melt(toReshape.fdi.px, id.vars = "x")
toReshape2.fdi.px <- data.frame(x = s.out.fdi.px[[2]][, "log(pfdidep_book)"], t(s.out.fdi.px[[6]]$fd) + t(s.out.fdi.px[[6]]$ev))
longFD.fdi.px <- melt(toReshape2.fdi.px, id.vars = "x")
longEV.fdi.px$Setting <- "Low Chinese FDI Dependence on Partner"
longFD.fdi.px$Setting <- "High Chinese FDI Dependence on Partner"
longEV.fdi.px <- data.frame(rbind(longEV.fdi.px, longFD.fdi.px))
longEV.fdi.px$x_percent <- exp(longEV.fdi.px$x)*100

# Plot
zp.fdi.px.facet <- ggplot(longEV.fdi.px, aes(x = x_percent, y = value, group = variable)) +
  geom_line(alpha = I(1/sqrt(nrow(s.out.fdi.px[[6]]$ev)))) +
  facet_grid(Setting~., scales="free_y", labeller = "label_value") +
  scale_x_log10("Partner's FDI Dependence on China (%), scale = log_10",
                breaks=c(0.1, 1, 5, 10), 
                labels=c("0.1", "1", "5", "10"), 
                expand = c(0, 0)) +
  scale_y_continuous("Predicted Probability of BSAs", limits = c(0, 1), expand = c(0, 0.01)) +
  theme_bw()
pdf(file="predictedprob_facet_fdi_px_bw.pdf", width=7, height=7)
print(zp.fdi.px.facet)
dev.off()

##############################################################################
## Fig 9. First Difference Estimates with Simulated 95% Confidence Intervals
##############################################################################
# combine simulated effect estimates
effect.df <- rbind(s.out.ptradedep.fd.stats, s.out.chtradedep.fd.stats, 
                   s.out.pfdidep_book.fd.stats, s.out.chfdidep_book.fd.stats, 
                   s.out.pta.fd.stats, s.out.bits.fd.stats, s.out.sco.fd.stats, 
                   s.out.oil_prod_eia.fd.stats, s.out.coal_prod_eia.fd.stats, 
                   s.out.gdp2000.fd.stats, s.out.rgdpch.fd.stats, 
                   s.out.gdpgrow.fd.stats, s.out.dist_k.fd.stats)
rownames(effect.df) <- c("ptradedep", "chtradedep", "pfdidep_book", "chfdidep_book",
                         "pta", "bits", "sco", "oil_prod_eia", "coal_prod_eia",
                         "gdp2000", "rgdpch", "gdpgrow", "dist_k") 
effect.vec <- effect.df[,1]
effect.low.vec <- effect.df[,3] 
effect.high.vec <- effect.df[,4]
effect.names <- var.names[1:13]

# prepare parameters for plot
y.axis.effect <- c(length(effect.vec):1)#create indicator for y.axis, descending so that R orders vars from top to bottom on y-axis

pdf("effectplot.pdf", height = 7, width = 7)#open pdf device
par(mfrow=c(1, 1))#reset graphical window
par(mar=c(2, 12, .5, 0))#set margins for plot, leaving lots of room on left-margin (2nd number in margin command) for variable names
plot(effect.vec, y.axis.effect, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = .6, col = "black",#plot coefficients as points, turning off axes and labels. 
     xlim = c(-0.9, 0.9), ylim = c(min(y.axis.effect), max(y.axis.effect)), xaxs = "r", main = "") #set limits of x-axis so that they include mins and maxs of 
#coefficients + 95% confidence intervals and plot is symmetric; use "internal axes", and leave plot title empty
segments(effect.low.vec, y.axis.effect, effect.high.vec, y.axis.effect, lwd =  1.5, col = "black")#coef +/-1.96*se = 95% interval, lwd adjusts line thickness
axis(1, at = seq(-0.9, 0.9, by = 0.1), labels =  seq(-0.9, 0.9, by = 0.1), tick = T,#draw x-axis and labels with tick marks
     cex.axis = .8, mgp = c(2, .5, 0))#reduce label size, moves labels closer to tick marks
axis(2, at = y.axis.effect, label = effect.names, las = 1, tick = T, 
     cex.axis = .8) #draw y-axis with tick marks, make labels perpendicular to axis and closer to axis
segments(0, 0, 0, 17, lty = 2) #draw dotted line through 0
dev.off()


################################################################################
## Table 9. First Difference Estimates with Simulated 95% Confidence Intervals
################################################################################
effects.table <- data.frame(effect.names, effect.vec, effect.low.vec, effect.high.vec)
names(effects.table) <- c("Covariate","Estimates","2.5%","97.5%")
rownames(effects.table) <- NULL
effects.xtable <- xtable(effects.table, caption = "First Difference Estimates with Simulated 95\\% Confidence Intervals", 
                         label = "table:first_difference", digits = 3)
print.xtable(effects.xtable, floating = TRUE, table.placement = "!ht", caption.placement = "bottom", 
             include.rownames = FALSE, include.colnames = TRUE, booktabs = TRUE)


#########
## PCPs
#########
## average fitted values for each model 
# de facto
fv.1 <- (relogit.robust.1[[1]]$fitted.values + relogit.robust.1[[2]]$fitted.values + relogit.robust.1[[3]]$fitted.values + 
           relogit.robust.1[[4]]$fitted.values + relogit.robust.1[[5]]$fitted.values + relogit.robust.1[[6]]$fitted.values + 
           relogit.robust.1[[7]]$fitted.values + relogit.robust.1[[8]]$fitted.values + relogit.robust.1[[9]]$fitted.values + 
           relogit.robust.1[[10]]$fitted.values)/10
# de jure
fv.2 <- (relogit.robust.2[[1]]$fitted.values + relogit.robust.2[[2]]$fitted.values + relogit.robust.2[[3]]$fitted.values + 
           relogit.robust.2[[4]]$fitted.values + relogit.robust.2[[5]]$fitted.values + relogit.robust.2[[6]]$fitted.values + 
           relogit.robust.2[[7]]$fitted.values + relogit.robust.2[[8]]$fitted.values + relogit.robust.2[[9]]$fitted.values + 
           relogit.robust.2[[10]]$fitted.values)/10
# full model
fv.7 <- (relogit.robust.7[[1]]$fitted.values + relogit.robust.7[[2]]$fitted.values + relogit.robust.7[[3]]$fitted.values + 
           relogit.robust.7[[4]]$fitted.values + relogit.robust.7[[5]]$fitted.values + relogit.robust.7[[6]]$fitted.values + 
           relogit.robust.7[[7]]$fitted.values + relogit.robust.7[[8]]$fitted.values + relogit.robust.7[[9]]$fitted.values + 
           relogit.robust.7[[10]]$fitted.values)/10

# generating the model predictions and a prediction success table.
fv.7.long <- t(fv.7)
predict.value <- ifelse(fv.7.long > 0.5, 1, 0)
data.predict <- merge(data, predict.value, by = c("row.names"), all = TRUE)

# rename variable
data.predict <- rename(data.predict, c(V1="predict"))

# table
predict.table <- table(predict.value,data$swap_lead)
colnames(predict.table) <- c("Actual 0", "Actual 1")
rownames(predict.table) <- c("Predict 0", "Predict 1")
predict.table

# calculate pcp
pcp <- (predict.table[1, 1] + predict.table[2, 2])/(predict.table[1, 1] + predict.table[1, 2] + predict.table[2, 1] + predict.table[2, 2])
# pcp among actual 0
pcp.0 <- predict.table[1, 1]/(predict.table[1, 1]+predict.table[2, 1])
# pcp among actual 1
pcp.1 <- predict.table[2, 2]/(predict.table[2, 2]+predict.table[1, 2])

###############################################################
## Fig 10. Rare Events Logit Model Fit: ROC Curve Comparisons
###############################################################
pdf(file = "ROCplots.pdf", width = 7, height = 4)
par(mfrow = c(1, 2), mai = c(0.8, 0.5, 0.6, 0.5))

# Full vs. defacto 
rocplot(relogit.robust.7[[1]]$y, relogit.robust.1[[1]]$y, fv.7, fv.1)
segments(0.1, 0.2, 0.15, 0.2, lwd = 1, col = "black")
text(0.15, 0.2, "Full Model", pos = 4, cex  = 0.8)
segments(0.1, 0.15, 0.15, 0.15, lwd = 1, lty = "dashed", col = "black")
text(0.15, 0.15, "De facto Model", pos = 4, cex  = 0.8)
# Full vs. dejure
rocplot(relogit.robust.7[[1]]$y, relogit.robust.2[[1]]$y, fv.7, fv.2)
segments(0.1, 0.2, 0.15, 0.2, lwd = 1, col = "black")
text(0.15, 0.2, "Full Model", pos = 4, cex  = 0.8)
segments(0.1, 0.15, 0.15, 0.15, lwd = 1, lty = "dashed", col = "black")
text(0.15, 0.15, "De jure Model", pos = 4, cex  = 0.8)
dev.off()

########################################################################################
## Robustness Check Tables
########################################################################################
# copper model
coef.names.copper <- rownames(coef(full.sum.copper))
coef.values.copper <- coef(full.sum.copper)[,1]
se.values.copper <- coef(full.sum.copper)[,2]
p.val.copper <- coef(full.sum.copper)[,4]
aic.copper <- (full.out.copper[[1]][[11]] + full.out.copper[[2]][[11]] + full.out.copper[[3]][[11]] +
                 full.out.copper[[4]][[11]] + full.out.copper[[5]][[11]] + full.out.copper[[6]][[11]] +
                 full.out.copper[[7]][[11]] + full.out.copper[[8]][[11]] + full.out.copper[[9]][[11]] +
                 full.out.copper[[10]][[11]])/10

# iron model
coef.names.iron <- rownames(coef(full.sum.iron))
coef.values.iron <- coef(full.sum.iron)[,1]
se.values.iron <- coef(full.sum.iron)[,2]
p.val.iron <- coef(full.sum.iron)[,4]
aic.iron <- (full.out.iron[[1]][[11]] + full.out.iron[[2]][[11]] + full.out.iron[[3]][[11]] + 
               full.out.iron[[4]][[11]] + full.out.iron[[5]][[11]] + full.out.iron[[6]][[11]] +
               full.out.iron[[7]][[11]] + full.out.iron[[8]][[11]] + full.out.iron[[9]][[11]] +
               full.out.iron[[10]][[11]])/10

# potash model
coef.names.potash <- rownames(coef(full.sum.potash))
coef.values.potash <- coef(full.sum.potash)[,1]
se.values.potash <- coef(full.sum.potash)[,2]
p.val.potash <- coef(full.sum.potash)[,4]
aic.potash <- (full.out.potash[[1]][[11]] + full.out.potash[[2]][[11]] + full.out.potash[[3]][[11]] +
                 full.out.potash[[4]][[11]] + full.out.potash[[5]][[11]] + full.out.potash[[6]][[11]] +
                 full.out.potash[[7]][[11]] + full.out.potash[[8]][[11]] + full.out.potash[[9]][[11]] +
                 full.out.potash[[10]][[11]])/10

# oilexp model
coef.names.oilexp <- rownames(coef(full.sum.oilexp))
coef.values.oilexp <- coef(full.sum.oilexp)[,1]
se.values.oilexp <- coef(full.sum.oilexp)[,2]
p.val.oilexp <- coef(full.sum.oilexp)[,4]
aic.oilexp <- (full.out.oilexp[[1]][[11]] + full.out.oilexp[[2]][[11]] + full.out.oilexp[[3]][[11]] + 
                 full.out.oilexp[[4]][[11]] + full.out.oilexp[[5]][[11]] + full.out.oilexp[[6]][[11]] +
                 full.out.oilexp[[7]][[11]] + full.out.oilexp[[8]][[11]] + full.out.oilexp[[9]][[11]] +
                 full.out.oilexp[[10]][[11]])/10

# coalexpwt model
coef.names.coalexp <- rownames(coef(full.sum.coalexp))
coef.values.coalexp <- coef(full.sum.coalexp)[,1]
se.values.coalexp <- coef(full.sum.coalexp)[,2]
p.val.coalexp <- coef(full.sum.coalexp)[,4]
aic.coalexp <- (full.out.coalexp[[1]][[11]] + full.out.coalexp[[2]][[11]] + full.out.coalexp[[3]][[11]] +
                  full.out.coalexp[[4]][[11]] + full.out.coalexp[[5]][[11]] + full.out.coalexp[[6]][[11]] +
                  full.out.coalexp[[7]][[11]] + full.out.coalexp[[8]][[11]] + full.out.coalexp[[9]][[11]] +
                  full.out.coalexp[[10]][[11]])/10

# exclude HKG
coef.names.noHKG <- rownames(coef(full.sum.noHKG))
coef.values.noHKG <- coef(full.sum.noHKG)[,1]
se.values.noHKG <- coef(full.sum.noHKG)[,2]
p.val.noHKG <- coef(full.sum.noHKG)[,4]
aic.noHKG <- (full.out.noHKG[[1]][[11]] + full.out.noHKG[[2]][[11]] + full.out.noHKG[[3]][[11]] + 
                full.out.noHKG[[4]][[11]] + full.out.noHKG[[5]][[11]] + full.out.noHKG[[6]][[11]] +
                full.out.noHKG[[7]][[11]] + full.out.noHKG[[8]][[11]] + full.out.noHKG[[9]][[11]] +
                full.out.noHKG[[10]][[11]])/10
n.noHKG <- length(mi1.noHKG$iso3)

# exclude MAC
coef.names.noMAC <- rownames(coef(full.sum.noMAC))
coef.values.noMAC <- coef(full.sum.noMAC)[,1]
se.values.noMAC <- coef(full.sum.noMAC)[,2]
p.val.noMAC <- coef(full.sum.noMAC)[,4]
aic.noMAC <- (full.out.noMAC[[1]][[11]] + full.out.noMAC[[2]][[11]] + full.out.noMAC[[3]][[11]] + 
                full.out.noMAC[[4]][[11]] + full.out.noMAC[[5]][[11]] + full.out.noMAC[[6]][[11]] +
                full.out.noMAC[[7]][[11]] + full.out.noMAC[[8]][[11]] + full.out.noMAC[[9]][[11]] +
                full.out.noMAC[[10]][[11]])/10
n.noMAC <- length(mi1.noMAC$iso3)

# export interdependence
coef.names.exdep <- rownames(coef(full.sum.exdep))
coef.values.exdep <- coef(full.sum.exdep)[,1]
se.values.exdep <- coef(full.sum.exdep)[,2]
p.val.exdep <- coef(full.sum.exdep)[,4]
aic.exdep <- (full.out.exdep[[1]][[11]] + full.out.exdep[[2]][[11]] + full.out.exdep[[3]][[11]] +
                full.out.exdep[[4]][[11]] + full.out.exdep[[5]][[11]] + full.out.exdep[[6]][[11]] +
                full.out.exdep[[7]][[11]] + full.out.exdep[[8]][[11]] + full.out.exdep[[9]][[11]] +
                full.out.exdep[[10]][[11]])/10

# import interdependence
coef.names.imdep <- rownames(coef(full.sum.imdep))
coef.values.imdep <- coef(full.sum.imdep)[,1]
se.values.imdep <- coef(full.sum.imdep)[,2]
p.val.imdep <- coef(full.sum.imdep)[,4]
aic.imdep <- (full.out.imdep[[1]][[11]] + full.out.imdep[[2]][[11]] + full.out.imdep[[3]][[11]] + 
                full.out.imdep[[4]][[11]] + full.out.imdep[[5]][[11]] + full.out.imdep[[6]][[11]] + 
                full.out.imdep[[7]][[11]] + full.out.imdep[[8]][[11]] + full.out.imdep[[9]][[11]] + 
                full.out.imdep[[10]][[11]])/10

# pimp.chexp
coef.names.pimp.chexp <- rownames(coef(full.sum.pimp.chexp))
coef.values.pimp.chexp <- coef(full.sum.pimp.chexp)[,1]
se.values.pimp.chexp <- coef(full.sum.pimp.chexp)[,2]
p.val.pimp.chexp <- coef(full.sum.pimp.chexp)[,4]
aic.pimp.chexp <- (full.out.pimp.chexp[[1]][[11]] + full.out.pimp.chexp[[2]][[11]] + full.out.pimp.chexp[[3]][[11]] +
                     full.out.pimp.chexp[[4]][[11]] + full.out.pimp.chexp[[5]][[11]] + full.out.pimp.chexp[[6]][[11]] +
                     full.out.pimp.chexp[[7]][[11]] + full.out.pimp.chexp[[8]][[11]] + full.out.pimp.chexp[[9]][[11]] +
                     full.out.pimp.chexp[[10]][[11]])/10

# pexp.chimp
coef.names.pexp.chimp <- rownames(coef(full.sum.pexp.chimp))
coef.values.pexp.chimp <- coef(full.sum.pexp.chimp)[,1]
se.values.pexp.chimp <- coef(full.sum.pexp.chimp)[,2]
p.val.pexp.chimp <- coef(full.sum.pexp.chimp)[,4]
aic.pexp.chimp <- (full.out.pexp.chimp[[1]][[11]] + full.out.pexp.chimp[[2]][[11]] + full.out.pexp.chimp[[3]][[11]] +
                     full.out.pexp.chimp[[4]][[11]] + full.out.pexp.chimp[[5]][[11]] + full.out.pexp.chimp[[6]][[11]] +
                     full.out.pexp.chimp[[7]][[11]] + full.out.pexp.chimp[[8]][[11]] + full.out.pexp.chimp[[9]][[11]] +
                     full.out.pexp.chimp[[10]][[11]])/10

# clean d.prob vectors
nam.d.prob.7 <- rownames(coef(d.prob.sum.7))
coef.d.prob.7 <- coef(d.prob.sum.7)[,1]
se.d.prob.7 <- coef(d.prob.sum.7)[,2]
pval.d.prob.7 <- coef(d.prob.sum.7)[,4]
remove.d.prob.7.table <- c(2, 5, 6, 9, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32)
nam.d.prob.7.table <- nam.d.prob.7[-remove.d.prob.7.table]
coef.d.prob.7.table <- coef.d.prob.7[-remove.d.prob.7.table]
se.d.prob.7.table <- se.d.prob.7[-remove.d.prob.7.table]
pval.d.prob.7.table <- pval.d.prob.7[-remove.d.prob.7.table]
aic.d.prob <- (d.prob.7[[1]][[11]] + d.prob.7[[2]][[11]] + d.prob.7[[3]][[11]] + d.prob.7[[4]][[11]] + 
                 d.prob.7[[5]][[11]] + d.prob.7[[6]][[11]] + d.prob.7[[7]][[11]] + d.prob.7[[8]][[11]] +
                 d.prob.7[[9]][[11]] + d.prob.7[[10]][[11]])/10

# prepare for texreg
full.copper <- MI(names=coef.names.copper, 
                  coef=coef.values.copper, 
                  se=se.values.copper, 
                  pval=p.val.copper, 
                  n=n,
                  aic=aic.copper)
full.iron <- MI(names=coef.names.iron, 
                coef=coef.values.iron, 
                se=se.values.iron, 
                pval=p.val.iron, 
                n=n,
                aic=aic.iron)
full.potash <- MI(names=coef.names.potash, 
                  coef=coef.values.potash, 
                  se=se.values.potash, 
                  pval=p.val.potash, 
                  n=n,
                  aic=aic.potash)
full.oilexp <- MI(names=coef.names.oilexp, 
                  coef=coef.values.oilexp, 
                  se=se.values.oilexp, 
                  pval=p.val.oilexp, 
                  n=n,
                  aic=aic.oilexp)
full.coalexp <- MI(names=coef.names.coalexp, 
                   coef=coef.values.coalexp, 
                   se=se.values.coalexp, 
                   pval=p.val.coalexp, 
                   n=n,
                   aic=aic.coalexp)
full.noHKG <- MI(names=coef.names.noHKG, 
                 coef=coef.values.noHKG, 
                 se=se.values.noHKG, 
                 pval=p.val.noHKG, 
                 n=n.noHKG,
                 aic=aic.noHKG)
full.noMAC <- MI(names=coef.names.noMAC, 
                 coef=coef.values.noMAC, 
                 se=se.values.noMAC, 
                 pval=p.val.noMAC, 
                 n=n.noMAC,
                 aic=aic.noMAC)
full.exdep <- MI(names=coef.names.exdep, 
                 coef=coef.values.exdep, 
                 se=se.values.exdep, 
                 pval=p.val.exdep, 
                 n=n,
                 aic=aic.exdep)
full.imdep <- MI(names=coef.names.imdep, 
                 coef=coef.values.imdep, 
                 se=se.values.imdep, 
                 pval=p.val.imdep, 
                 n=n,
                 aic=aic.imdep)
full.pimp.chexp <- MI(names=coef.names.pimp.chexp, 
                      coef=coef.values.pimp.chexp, 
                      se=se.values.pimp.chexp, 
                      pval=p.val.pimp.chexp, 
                      n=n,
                      aic=aic.pimp.chexp)
full.pexp.chimp <- MI(names=coef.names.pexp.chimp, 
                      coef=coef.values.pexp.chimp, 
                      se=se.values.pexp.chimp, 
                      pval=p.val.pexp.chimp, 
                      n=n,
                      aic=aic.pexp.chimp)
d.prob <- MI(names=nam.d.prob.7.table, 
             coef=coef.d.prob.7.table, 
             se=se.d.prob.7.table, 
             pval=pval.d.prob.7.table, 
             n=n,
             aic=aic.d.prob)
# Table 7. Robustness Check Models I: Parameter Estimates and Robust Standard Errors
texreg(list(full.oilexp, full.coalexp, full.copper, full.iron, full.potash, full.noHKG, full.noMAC),
       stars = c(0.001, 0.01, 0.05, 0.1), symbol = "\\wedge", sideways = TRUE, single.row = TRUE, 
       custom.model.names = c("Oil Exp.","Coal Exp.","Copper","Iron","Potash","Exclude HKG","Exclude MAC"), 
       custom.coef.names = c("Intercept", "log(P. Trade Dependence)","log(CN Trade Dependence)","log(P. FDI Dep.)",
                             "log(CN FDI Dep.)","PTAs","BITs","SCO","log(Oil Exports)","log(Coal Production)",
                             "log(GDP)","log(GDP per capita)","GDP Growth Rate","Bilateral Distance",
                             "log(P. Trade Dep.):log(CN Trade Dep.)","log(P. FDI Dep.):log(CN FDI Dep.)",
                             "log(Oil Production)","log(Coal Exports)","log(Copper Exports)","log(Iron Exports)",
                             "log(Potash Exports)"), 
       reorder.coef= c(1, 2, 3, 15, 4, 5, 16, 6, 7, 8, 17, 9, 10, 18, 11, 12, 13, 14, 19, 20, 21), float.pos = "!ht", 
       caption = "Robustness Check Models I: Parameter Estimates and Robust Standard Errors", 
       label = "table:robust_resources_sar", use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, scriptsize = TRUE)

# Table 8. Robustness Check Models II: Parameter Estimates and Robust Standard Errors
texreg(list(full.exdep, full.imdep, full.pimp.chexp, full.pexp.chimp, d.prob), 
       stars = c(0.001, 0.01, 0.05, 0.1), symbol ="\\wedge", sideways = TRUE, single.row=TRUE, 
       custom.model.names = c("Exp. Interdep.","Imp. Interdep.","P. Imp*CN Exp","P. Exp.*CN Imp.","D.Probit"), 
       custom.coef.names = c("Intercept","log(P. Exp. Dep.)","log(CN Exp. Dep.)","log(P. Import Dep.)",
                             "log(CN Import Dep.)","log(P. FDI Dep.)","log(CN FDI Dep.)","PTAs","BITs","SCO",
                             "log(Oil Production)","log(Coal Production)","log(GDP)","log(GDP per capita)",
                             "GDP Growth Rate","Bilateral Distance","log(P. Exp. Dep.):log(CN Exp. Dep.)",
                             "log(P. FDI Dep.):log(CN FDI Dep.)","log(P. Imp. Dep.):log(CN Imp. Dep.)",
                             "log(P. Imp. Dep.):log(CN Exp. Dep.)","log(P. Exp. Dep.):log(CN Imp. Dep.)",
                             "log(P. Trade Dependence)","log(CN Trade Dependence)","log(P. Trade Dep.):log(CN Trade Dep.)"), 
       reorder.coef= c(1, 22, 23, 24, 2, 3, 4, 5, 17, 19, 20, 21, 6, 7, 18, 8, 9, 10, 11, 12, 13, 14, 15, 16), float.pos = "!ht", 
       caption = "Robustness Check Models II: Parameter Estimates and Robust Standard Errors", 
       label = "table:robust_dis_trade", use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE, scriptsize = TRUE)

